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ABSTKACT 


AUTOCORRELATION AND POWER SPECTRUM ANALYSIS TOR 
X-RAY AND GAMMA RAY SPECTROMETER DATA 

by 

CHUG- JEN TSAI 

The purpose of this thesis is to test the usefulness of autocorre- 
lation and power spectrum analysis computer programs for studying signals 
from possible X-ray and gamma ray pulsar emitting pulsars such as NP 0532. 
For checking the program, simulated data are generated: (1) simulated 

random background, (2) simulated square pulse signal above a random back- 
ground and (3) a simulated pulsar signal above a random background. These 
data are discrete and equally spaced time series. 

The results of the analysis, when the simulated pulsar signal repre- 
sents a 100% intensity increase above the random background, is in good 
agreement with an analytical solution. This technique fails to detect pul- 
sar signal if it has an intensity less than 25% above the random background 
_ 5 

for a set of 10 data points. The sensitivity of the technique will be 
improved for an increased amount of data or increased observation time. 

Also if the period is known the parameters used in the analysis may be opti 
mized to increase the sensitivity of the method. 


ix 



CHAPTER I 


INTRODUCTION' ' 

1.1 General Aim of Thesis 

The emission of X-rays from the pulsar NP 0532 in the Crab Nebula 
has been established by the Naval Research Laboratory Group (FRIEDMAN et 
al. , 1969) by using the autocorrelation and power spectrum analysis tech- 
nique. Similarly, X— ray or gamma ray pulsations of NP 0532 could be ob- 
served by power spectrum analysis for balloon data, where flight times 
are several hours long and the random noise background is much larger than 
in rocket flights because of the atmospheric photons. 

The purpose of this thesis is to test the usefulness of a computer 
program to calculate the autocorrelation function and power spectrum with 
a simulated pulsar signal superimposed' on a random background. Before 
testing the program for calculating the autocorrelation function and power 
spectrum, three kinds of simulated data are generated. These are: 

(1) simulated random background, (2) a simulated square pulse signal on 
the random background, and (3) a simulated pulsar signal on the random 
background. The first two are subsidiary data for checking the autocorre- 
lation function and power spectrum computer programs. 
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1.2 Review of Observations of Pulsar NP 0532 


Pulsars, radio-emitting "stars'* having rapid and extremely accu- 
rate repetitive changes in luminosity with time, were discovered in 1968 
by HEWISH, et al. (1968). Later, STAELIN and KEIFENSTEIN (1968) observed 
two radio pulsars in the vicinity of the Crab Nebula. Of the-, pulsars yet 
discovered, these two have the longest (3.745 s.) and shortest periods 
(33.09 ms.). 

In January of 1969, COCKE, DISNEY, and TAYLOR (1969) reported the 
discovery of optical light flashes from the Crab Nebula. The flashes 
occurred with the same periodicity as the fast Crab pulsar (NP 0532) and 
were suggested to have originated south of the two central stars in the 
Crab Nebula. 

X-ray pulsations from NP 0532 as reported by NEL (FRIEDMAN et al., 
1969) indicate that it pulsates at a frequency closely matching the radio 
and optical pulsations. About 5% of the total average X-ray power of the 
Nebula appears in the pulsed component. The X-ray pulsations have the 
form of a main pulse and an inter-pulse, separated by about 12 ms. The 
important characteristics of NP 0532 (FRIEDMAN et al., 1969) for radio, 
optical and X-ray observations are tabulated in Table 1. 
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Table 1 


Radio Optical X-ray 


Average Pulsed 
Power 

(erg cm ^ sec 

6 x 10~ 14 
(195-430 Mhz) 

• -12 
8 x 10 

(4500-8500A°) 

1.5 x 10“ 9 
(1-10A°) 

Spectral 

Index 

-2 

• * 

U-B = -1.3 
B-V = +0.1 

0.4 

Separation of 
Main Pulse 
and Inter pulse 

'v 14.5 ms 

14.0 ms 

12 .0 ms 

" 

Half-power Width 
of Main Pulse 

^ 3.0 nis 

1.4 ms 

i 2. 5 ms 

Half-power Width 
of Interpulse 


3.0 ms 

< _ - 

= 5. 0 ms 


* U, B, and V are logarithmic intensities (magnitudes) in the ultravio- " 
let, blue, and visual, respectively. 
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1.3 The Reasons fo.r Using Autocorrelation and Power Spectrum Analysis 
' for Determining .the Pulsar Periods 

Autocorrelation and power spectrum analysis reveals information on 
pulsar characteristics in the time and frequency ‘domain. Suppose that an 
X-ray or gamma ray pulsar has one pulsation period, and the intensity of 
the signal is strong with respect to the random background, then the auto- 
correlation function will be a periodic function with the same period as 
the signal. In the frequency domain, the power spectrum would have a 
large value at several definite frequencies corresponding to the pulsar’s 
fundamental frequency, the first harmonic frequency, etc. 

In general, there might he several frequencies in the pulsar in 
question or in other words, it might have several pulsation periods. In 
this situation the autocorrelation function in the lag time domain is not 
a good way to pick up the pulsation periods. Because the autocorrelation 
can be considered as the sum of the- autocorrelation functions for the 
separate periods superimposed on the random background, the resultant auto- 
correlation function in the lag time domain might not be easily interpre- 
ted. On the contrary, the power spectrum allows one to easily detect those 
frequency components present in the pulsar. So long as this spectrum, when 
plotted versus frequency, has bumps or large values at several particular 
frequencies, one can easily tell which frequencies are the fundamental pul— 

■ sation frequencies. 

The power spectrum analysis determines not only the estimated value 
of the period but, it can determine how much power is contributed by each 
frequency component. The power contributed by each component is proportional 
to the square of the amplitude of that component. 
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1.4 Statement of Problem . 

Knowing the actual characteristics of a typical pulsar such as 
NP 0532, we could simulate this data for checking the "autocorrelation 
and power spectrum" program. Using instead an idealized pulsar simula- 
tion we can study several questions relating to the use of autocorrela- 
tion and power spectrum analysis techniques. The basic questions studied 
in this thesis are: (1) what is the limitation of using the autocorrela- 
tion and power spectrum analysis techniques as the pulsar signal becomes 
weak compared to the background? (2) What are the various advantages 
and disadvantages of using the technique? 
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CHAPTER II 
THEORY 

2.1 Theory of Autocorrelation and Power Spectrum Analysis 

Autocorrelation and power spectrum analysis is a very useful tool 
for detecting periodic signals buried in noise and for establishing coher- 
ence between random signals. Its applications range from engineering to 
radar and astronomy to medical, nuclear, and acoustical research. 

(1) Autocorrelation Function 

Autocorrelation for any kind of waveform is a measure of the rein- 
forcement between two identical waveforms shifted in phase. It Is compu- 
ted by multiplying one waveform ordinate by ordinate with the other and 
finding the average value of the product. If there is no phase shift 
involved, the correlation between two identical waveforms is large. In 
other words, the autocorrelation function at the zero lag time is large 
and a maximum. However, if two waveforms are identical in shape but have 
an arbitrary time shift between- them, then the correlation between them Is 
generally small. Hence, autocorrelation is a function of the time shift 
between the two identical waveforms. This correlation problem is illustra- 
ted in Figures 1 and 2. Figure 1 shows two identical waveforms without any 
time lag or phase shift. Figure 2 illustrates two identical waveforms with 
a lag time or phase shift. 

This function can be expressed in terms of a mathematical formula as 
the product of the wave x(t) and a delayed version of itself x(t + t) aver- 




• FIGURE 1 
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aged over T seconds as follows: 

= Um J- f'xctj xct + r> H ( 1 ) 

co l o 

where x is the lag time. This quantity R(t) is always a real-valued even 
function with a maximum at x=0 and may be either positive, negative or 
zero for other values of x. There may be other maxima depending on the 
function. Therefore, this function has: 

1. Symmetry about x=0, i.e. R(x) = R(-x) (2) 

2. A maximum at x=0 equal to the mean square value (X 2 ) (3) 

of the signal from which it is derived, i.e. R(0) = X 2 

and R(0) > R(x) for all x. 

Also for the special case of a periodic waveform the autocorrelation func- 
tion is periodic and has the same period as the waveform itself. This peri- 
odic autocorrelation has the two properties listed above but it will have 
maxima whenever the lag time x is equal to an integral number of periods of 
the signal. 

The random noise signal is quite different from the periodic wave- 
form. When compared with a time shifted version of itself, only a small 
time shift is required to destroy the correlation, and it never recurs. 

The autocorrelation function for this case is, therefore, a sharp impulse 
which decays from the central maximum to low values at large time shifts. 

Two samples of random noise of the same bandwidth might have quite 
different waveforms, but their autocorrelation functions could be identi- 
cal. The autocorrelation function of any signal, random or periodic, 
depends not on the actual waveform but on its frequency content. 
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(2) Power Spectrum 

The power spectral density function at frequency f £ is defined by 


BENDAT (1958) as: 


*<**.*) 'g. m ^ > 


( 4 ) 


where P(f c> x » Af) is the total average power in a given bandwidth Af and 


x is a continuous variable. Therefore, the above equation represents the 
limit of the total average power in a given bandwidth divided by the band- 
width as the bandwidth approaches zero. Suppose x(t) represents an infi- 
nite record length, then we can define data with a finite record length as 

x(t) l.tl < T 


x T (t) - 


0 


otherwise 


The total average power of bandwidth Af is defined by BENDAT (1958) as: 

•A -f* 

" * IAr<tyl 


1 * T~)oo J. 


JC. JL. 


~T 




(5) 


where 


h Of, %) = ( z(Ve. dt=j z -r a > e tit 

1 -T 


~<x> 


(6) 


^(fjx) is the direct Fourier transformation of the finite length of data 
^(t) and 

|Ar<f,w|“ A*tf.v A T (i, tj ■ (7) 

From equations (4) and (5) we observe that since |A,j(f,x)| 2 is an even 
function of f, it follows that G(f,x) is an even function for all f. The 
notation G(f,x) shows clearly that the power spectral density function de- 
pends upon the particular time record x(t) under examination. 
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2.2 The Relation Between the Autocorrelation Function 


and the Power Spectrum 


Consider an arbitrary real-valued time record x(t) of infinite 

extent. Its total energy is defined by 

oi 


r -a. 

Total energy = I K ^ ) ^~t 


( 8 ) 


— oQ 


If x(t) does not approach zero rapidly enough for large values of t, the 
total energy may be infinite; that is, this integral fails to converge. 
However, we shall assume that the average power associated with x(t) is 
finite, where the average power is defined by 

T Ai> ^ < s > 


T~> 


CO 


-T 


Using Parseval's Theorem as given by KHARKEVICH (1960a). let A^(f,x^) 
and A 2 (f,X 2 ) be the Fourier transformation of real functions x^(t) and 
X 2 (t) , respectively. Then 

J %{-t)%&) dt = J At Z >) fit X <*> (1Q) 




For the special case of x^(t) = ^(t) ain ^ a record 

f Xjitjat =J j A t a, %) f V = .2 J IA. T (*,xjfa 

— rA ^ 


If 


( 

Xj(-b)Ctt is convergent, then by (11) 

. r = L! „ 


^ i l 


T 




T 


if 


(ii) 


( 12 ) 



12 


Therefore, from equations (3), (4), (5), (11) and (12) 

fcr O) ~ xh-t) = Urn r x\*)cl± = f\dJ 4f 

”7" j* J 

I C> 


(13) 


but 


£(o;=J fc('?jSCc)d? 

-ob * 


because R(t) and 5(r) are even functions. 

Therefore, by use of the definition of the Dirac delta function <5(r) 
r °d r°* 

2 J d (TJ 4* 

O tL od. 


R(v 


dtj- 4' 


(14) 


Comparing (13) and (14) 


zrfr 


r oS —A* 1 J <■ 

2} R(-V£ * 


-ob 

Therefore, the power spectrum is the Fourier transform of the autocorre- 
lation function. The power spectrum can be written in another form as 
follows 


Rnjc°s **£?<*■ 


Since R(t) is an even function of t, G(f) is always a real-valued non- 
negative function. In short, the autocorrelation function and power spec- 
trim are the Fourier transformation pair, i.e., 

ftm) ~ p <JC fjc*s 

0 (15) 

Crtfj = 4 J ftcyeosinfrtf'* 

c 

Autocorrelation Theorem 

There is a useful theorem which is stated by BRACEWELL (1965) and 
can be used to explain the relation between power spectrum and autocorre- 
lation function. If X(t) has the Fourier transform A(f), then its auto- 
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correlation function is the Fourier transform of jA(f)| 2 . i.e.', if 

~ c zttJ t . . . 

j XCt)SL, dt .= 

K-«C> 


then. 


r 

— cO 


cri 


2 /' 2.7T-f'( J A~~ if- 

I Mf)j c. =J X tt) x.i-t + T)4t 




( 16 ) 


This autocorrelation function is unnormalized with zero mean. Thus 
I A (f ) | 2 is the power .spectrum in equation (16), since the Fourier trans- 


formation of the autocorrelation function is the power spectrum, that is 

Ri~)e- dn = $(+) 


i 


Conversely, the autocorrelation function is the Fourier transform of G(f) 


04 C^ir • 

J 6fih dj =RCrj K&) + 

-~c$ 


dt 


from the above equation compared with equation (16) we see that |A(f)| 2 


is equal to G(f). 
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2.3 Derivation of Autocorrelation Function and Power Spectrum 

for Computer Programming Use 


The functions used in the first two sections of this chapter were 
continuous for the convenience of the theoretical treatment. In practice 
digital type data are normally used. For example, the data obtained by 
collecting counts over some time interval At Is discrete. Therefore, in 
this sense, the simulated data we must use should be discrete. 

The simulated data will be equi-spaced, discrete and for a finite 
length record. For this reason, the formulas of autocorrelation function 
and power spectrum must be transformed. 

Suppose this series of data has a non-zero mean. An estimated 
autocorrelation function [I.B.M. Programmer’s Manual: System/360 

Scientific Subroutine Package, (360A -r CM 03X) , Version III, 59.] may be 
expressed as 


A 

R, * 


/ ... _ , _ 
- — - x i ) 


n~3 + l *=* 


(17) 


where: 


X is the mean value of all the data X., i.e. X = 

i 

n = number of observations in the time series X_^ 

j = 1, 2, 3, m represents time lags 0, 1, 2, (m-1) and 

m is the maximum lag number and the maximum lag time T - (m-1) x sampling 
time. For programming convenience, we may expand equation (17) into three 
terms as follows: 

yrjH n'j*/ 

Z 


I 


^ n-)*l 


r 


2 X. , - A Z LX- -t X. . ) +Z LX} 

C-i > '+J-I c*t J ~t ' ~C^I J 


(18) 
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The middle term can be expressed differently by using X according to the 
mathematical steps following: 


_ , >5 , n \ 

V r-LjX, = -(! -t 3 T X.) 

*1 c=S ° n •'-f * / 


• 

h'Jfl 


r) 


X 

=■ ” * - 

.ST. *, 


6-1 


*• = » 


x-Jtt 


n 

and 

T~ 

i~l 

X. . , = - 

rX K 
*0 K 



n 

, 5-/ 

but 

X : 

r ± r x. = 

n Cs/ c 

: -L/ X 


n 

2: 


V»*>! 




3-/ 


£ *£+*-, = -F \- 


t-l 


t=- 




(19) 


(20) 


Erom equation (19) and equation (20) 


n*3+/ 

XT 

<>/ 


.d'i - * i 

fr^.+ 2 ; 




( 21 ) 


Substituting equation (21) in equation (18) gives' the autocorrelation 
function with non-zero mean as follows 

? W-3+1 ( e-j ‘ "+ 3 -* U **-/ t-M-j+a 

~-*-(y>~j+l )x\ 


I 


»-j+i 


-- /?VT» J'/ I ,~t 

£ 5 ^ + * *J (22) 


K-j H 


Equation (22) is the complete expression for calculating the autocorrela- 
tion function with a non-zero mean for equi-spaced, discrete data. The 
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parenthesis of the second term of this equation is the total sum of the 
first (j-1) data samples and the sum of the last (j-l) data samples. 

BLACKMAN AND TUKEY (1958a), states that the ensemble average of the 
estimated autocorrelation function for the discrete case equals the true 
autocorrelation function for the continuous case. So that 


ave {Rj }■ = R(t) = R(jAx) 


Similarly, the average of the raw power spectral density can be expressed 


in the form of a Fourier transform, viz., 


} = 2 [ v„ Crj , rctjJ e rft 


where is a finite Dirac comb, viz.. 


j-m-1 




v (Q- - At) = -i- 4?r & fcnr 

w w w ~ * o=-«n 

Then = 2 J ±1 RC?) -h 




- L.«=h 

0 = 'hi-f| ■ — 

LtS^rilrdl: w . f • 

- ^ ^ CoS i/r-f + 3/e.J &>] 2T 

G>S2.n$ 


If we ignore the averaging (Ave) in the above equation, the raw spectral 
density may be expressed as: 

O'" A .A M-l <V 

Sr =. R„ <** x-foc** (23) 

n ° o’-/ J 

where f is not a continuous variable any more but is now a discrete par- 


ticular number. 



Let 


( J A t 

W^ere. 


o Kn 

ry\ 


e. 


*> ^ 


Then, equation (23) becomes, 


£ = -£— 
* £■ !» * ^ 


€p =i2.*~f R c + 2j 
A. , - ) ‘" 


^ jKn 

r 3 ^ ^ 



( 24 ) 


(25) 


In order to get a smoothed power spectrum the lag window and the 
spectral window are involved. The lag window is the Fourier transform of 
the spectral window. There are several Fourier transformation pairs 
(BLACKMAN and TUKEY, 1958b) relating these two windows: in the following 

D denotes lag window and Q denotes spectral window. 


Zeroth pair 


%(V 


= 1 


in < 7 - 


>M 


I'TI >Z 


and 


S/’w 


* 


2n f T „ 


First pair (BARTLETT, 1950) 

= i- V 


t~l < i 


v,m 




= O 


and 


Gli&J — Xh (- ~7rfr M ) 
Second pair (sometimes called "Hanning") 




>K 


- -‘-Cl-*- CoS ) 






1*1 < u 


=- o 


and 


1^1 >t m 

=f<a.aj + ^r Cfiot* +zfj -*•<?. rf ~Xrj] 

Third pair ( sometimes called "Hamming") 


D.CZ) 


-ss o.S‘4 *+ o-44 os 


— o 


7T < 


1*1 


m. 


and 


;~j > 7 ^ 
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Fourth pair 


\>4z. -t O' fro ~~ ~fr O'0$ 


2 >«^J - 


I'J'I 


and Q CT; = 0.42 a.(fi +o^C 4 «- r ^+t 3 .fJ L -£i^ 

f 4**4 £ So _Cf + ir ^ -f- <2. * yOj 

BENDAT and PIERSOL (1966a) use Hanning (second pair) for smoothing the 
raw power spectrum. 

The smoothed power spectrum is the convolution of Q (f) and the 

true power spectrum. Therefore 

X-\- 'v 

**- O'fr ^ -f- o.i - Cf t 

= ^ + ****** i **’'* ' (26) 

T = o,S 6r + o.fr 6f 

m-t 

KHARKEVICH (1960b) states that, "Any function f (t) consisting only 

of frequencies from 0 to f can, with any desired accuracy be treated as a 

succession of numbers recurring every ~~~ seconds." Where f .in the no- 

c 

tation of BLACKMAN and TUKEY is f , This maximum ‘frequency is also known 
as the folding (or Nyquist) frequency. From the equation (24), when K 
equals m then f ^ which is the Nyquist frequency. At in the notation 

of BENDAT and PIERSOL (1966b) is h or the time interval between samples. 

The resolution bandwidth for power spectrum is defined as B = ^ . The 

B will be small for a given h when m is large. 

An important feature known as "aliasing" enters for observing an 
equally spaced, discrete finite record of data. The energy or power at an 
arbitrary frequency f cannot in general be separated from that contributed 
by different frequencies [BLACKMAN and TUKEY (1958c)]. In other words, 
higher frequencies from the original process G(f) may contribute some power 
to the estimated power spectrum G .(f) (see Fig. 3). Figures 3a and 3b in- 

A 
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G(f.) 



'(a) 

G (f) 



FIGURE 3 
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dicate that the powers contributed by frequencies 2f -f^ and are indis- 
tinguishable. The essential, unavoidable nature of this problem is made 
clear by Fig. 4 which illustrates how equally spaced time samples from any 
cosine wave could have come from each of many other cosine waves. In 

Fig. 4, the sampling time Ax is 0.2 second; then the Nyquist frequency is 

1 1 

~ 2 ^ ~ "q—^ cycle/ sec . , i.e., 2.5 cycle/sec. Clearly, a sinusoidal wave with 
frequency 4 cycle/sec. can be fit through the points and a second sinusoi- 
dal wave (dotted curve) with a frequency 1 cycle/sec may also fit the given 
points. Higher frequencies may also be present but it is not possible to 
know from the measured At whether the power at frequency f [of a power spec- 
trum in the interval (0,f c )] comes from the principal frequency for 
2f -f, sf +f, 4f -f, 4f +f,.....etc. These higher frequencies are called 
aliases of f. Therefore, the aliased power spectrum G^(f) defined in the 
interval (0,f ) by BLACKMAN and TUKEY (1958d) is all that one may estimate 
from the data. 

G^(f) may be represented as follows: 

G A (f) - G(f ) + G(2f c -f) + G(2f c +f) + ...etc. 0< | f J <f ^ 

G A (f) = 0 

where G(f) is the true power spectrum. 

Two practical methods exist for handling this aliasing problem 
[BENDAT and PIERSOL (1966c) 3 • The first method is to choose h suffici- 
ently small so that it is physically unreasonable for data to exist above 
the associated Nyquist frequency f £ . For the low frequency case of inter- 
est, say below 500 HZ., then h = 1 ms would technically be sufficient. 

The second method is to filter the data prior to sampling so that no infor- 
mation above the Nyquist frequency is contained in the filtered data. Then 
choosing f as the maximum frequency of interest will give accurate results 
for frequencies below f c . 
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F.IGURE 4 
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CHAPTER III 
PROGRAMMING* 

The X-ray or Gamma ray spectrum of a pulsar might be in the form 
of a small signal on the random background. Autocorrelation and power 
spectrum analysis has already been used as a tool in conjunction with a 
rocket experiment (NRL Group) for detecting the X-ray pulsar (NP 0532) 
signal over the general Crab X-ray signal. Similarly, the autocorrelation 
and power spectrum technique might be used for balloon experiments. In 
order to check the programs for computing the autocorrelation and power 
spectral density function before applying them to balloon flight data, 
three kinds of simulated data are generated by the computer. The first 
one is a simulated random background for discrete data. The second one 
is a simulated square pulse signal above a random background. The third 
one is a simulated pulsar signal above the random background. 
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3.1 The Basic Principle Used for Generating 'Random Data 


A way of generating random data is to use the time interval dis- 
tribution which is derived from the Poisson distribution which gives the 
probability distribution in time for a random process when the average 
rate is specified. 

The time interval distribution describes the distribution in size 
of the time intervals between successive random counts when the mean rate 
has the constant value of a counts per unit time. From EVANS (1955) sup- 
pose a is the average rate of appearance of photons; then the average 
number of counts in a time interval t is at. 

' When the average rate is a, the probability of observing x counts 
in a time interval t can be expressed as follows: 


P 00 - e~ it; 

X - X. 


This is simply the Poisson Distribution. From this equation the 
probability that there will be no counts in a time interval t, during which 
time there should be at counts on the average, is 


p n = (at) ° e ~ at = e“ at 

u 0! 


and its differential form is shown as follows: 

dP 
o 

dt 

~We at once see that small time intervals between successive random 
counts have a higher probability (probability of no counts in a time inter- 
nal t) of occurring than large time intervals. In other words, this is the 
TIME INTERVAL DISTRIBUTION which gives the probability of occurrence of 



each time interval. 
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Of course, the occurrence of a given time interval is random and 
may be represented by a random number. This implies that the differential 
probability of occurrence of a given time interval may also be represented 
by a random number. 

Now from the random time interval between two successive counts, 
the random counts within a definite time interval could be generated. 

Those random counts are the simulated data I have been using. Suppose all 
events (counts) occur randomly along a time axis, the number of random 
counts is then nothing hut the number of counts occurring within a defi- 
nite time interval. 

Thus, the steps of generating random counts is shown as follows: 

1. Generate random number between 0 and 1 by use of the SUBROUTINE T RANDU. * 
(see Appendix A) 

2. Take this random number as the probability of no counts occurring in 

a time interval t. but one count in t. to t. + dt, and generate random 
x xxx 

r. time intervals according to the-f ormula t . = - . w ^ ere ^ is a 

a 

random time interval and UNIT is random number between 0 and 1. 

3. Generate the number of random counts within a given desired sampling 

time by observing the sequence of intervals obtained from step 2. The 

distribution of the random time intervals is shown in Figure 5 by using 

some of the data from the "Random Data" program. The random time inter- • 

val data have been generated under the average counting rate a=l count/ms. 

5 4 

The total number of random time intervals is 1.5x10 , but only 1.5x10 
are printed out with the form of 1500 numbers of random time intervals 
for each 10th record. In' other words, there are 1500' numbers of ran- 
dom time intervals for 10th, 20th, 30th up to 100th record. 

In order to show the distribution of time intervals, the grouping 
of the numbers of time 'interval in the range of 0 ^ 0.1 ms, , 0.1 ^ 0.2 ms, • 
etc. is tabulated in Table 2. 



Hange of Time 
Interval (ms) 

Numbers of Time 
Interval (counts) 

0.1 * 0.2 

150 

0.2 'v 0.3 

110 

0.3 n. 0.4 . 

102 

0.4 ^ 0.5 

84 

0.5 'v. 0.6 

82 

0.6 a 0.7 

72 

0.7 o, 0.8 

62 

0.8 0.9 

46 

0.9 n. 1.0 

62 

1.0 n, 1.1 

50 

1.1 * 1.2 

59 

1.2 Vl.3 

‘ 41 

1.3 ^ 1.4 

38 

1.4 -v 1.5 

33 

1.5 n, 1.6 

37 

1.6 m 1.7 

27 

1.7 'V 1.8 

40 

1.8 * 1.9 

24 

1.9 'v 2.0 

17 

2.0 'u 2.1 

21 


Total 

Humber 

(Counts) 


1173 


2.1 * 2.2 


16 
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There are only 1173 out of 1500 values of time intervals taken 
from the first 1500 values of printed data from the program of •'Genera- 
ting Random Data." 

Figure 5 is a semilogarithmic plot of the numbers of time inter- 
vals in a given sampling time of 1 ms versus the time interval. The 
straight line shows that the distribution is exponential as we expect. 

From the slope of the line, we can determine a for checking the "Random 
Data" program. 

As I mentioned above a = 1 count/ms gives 1 count within 1 ms on 
the average which is a very low value statistically. Therefore the distri- 
bution of x counts occur within 1 ms time interval will follow the Poisson 

C3- 1 

distribution ^ ; ■£• . Substituting this value in this 

equation, then ^(ij — —j ~ j- -4L ^ 5 where P (1) is a random number be- 

cause of x random counts. The slope of the straight line in Fig. 5 is 
1.17 cts/ms . The 17% error is due to the limited data (1173), this would 
decrease with more data. A least square fit may 'give a better result. . 
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3.2 Brief Statement of the Method for Generating a Square Pulse 
or a Pulsar Signal Including a Random Background 

From Section 3.1, we know that the random background has been gen- 
erated by three steps. ..Once the random background has been generated, the 
random counts for either the square pulse. or pulsar signal also could be 
generated by the same steps. 

The process of constructing simulated data can be understood by 
reference to Figure 6. This shows the case where in region I the random 
counts correspond to an average rate a and are generated by the method 
described in Section 3.1. In region II, the random counts are generated 
for a new average rate f times the former rate, where f is any positive 
real number. To construct the data for a continuous flow of time, the 
random counts for a are used for a time t^, and then for a time t ^ the ran- 
dom counts for the new average, fa, are used, then the counts for a again 
for t^ are used etc. This procedure constructs the time series for a 
pulse which is on for a time t^ and off for a time t^ so the period is 

t l + t 2* 

If we set t^ = t^ = 20 ms, and a = lct/ms, and fa = 2cts/ms, then 

the data would be the square pulse with a 100% intensity increase over the 

random background with period 40 ms. Similarly, a pulsar signal with a 
period of 40 ms and a 100% intensity above the random background could be 

generated in the same way by setting t^ = 35 ms, - 5 ms, a - lct/ms. 


and fa = 2cts/ms. 
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3.3 Flow Ghart for Autocorrelation. Function and 


Power Spectrum Computation 


We will now discuss the procedure for developing programs to cal- 
culate the autocorrelation function and power spectrum. It is most conven- 
ient to describe this in terms of a flow chart. 

The autocorrelation function with nonzero mean for discrete data is 

given by 

I n 'Jt ! / „ - * T. s 

(17) 


/?. - • — L — r tz c -x)t% -x) 


where X is the mean value of the data x., 

l 

X = — I x. 
n i=l 1 


n = number of observation in time series X. 

x 


3 = 1, 2, 3, 


m represents time lags 0, 1, 2, 


(m-1) 


The transformed equation 


A f 


o-t n _ , 

% K- *;.H J (22) 


The raw power spectral density function 

A A / 7TjlC , 


% c* / - % C & )= * k l ~ (^T) ^ K L} 


where f = 


Kfc 

m 


K = 0,. 1; 2 


m 


h is the sampling time interval 

A 

R, is the estimate of the autocorrelation- function at time lag j-1 
3 

av is the maximum lag numb er 

f = — 
c 2h 


is the Nyquist frequency 
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The smoothed power spectrum is given by 
= ©.££ i-O'S ^ 

§■ = % + O'** j K=> > 2 , •■• *-l 


(26) 


'k '«■ 

$ ~ otdf 

It is very cumbersome .to calculate the autocorrelation function throughout 


all the data. However, we can calculate it for one part of data and then 
go through the other part of data. Therefore, it is necessary to mention 
what the "read" process is. This process is divided into three steps as 
may be seen by referring to the flow chart shown in Figure 7. 

(1) Initial correction: Read first 500 data values and do calcu- - 

lations of 

Jroo , . 500 

X J iron I -to too ; /\\JZR= T >v 

(see Figure 7, blocks 12, 13, 14, 15, 16, 17, and 18) 

(2) Main loop: Read 1000 data values (i.e. half of the data in 

first record and half of the data in second record) and do the calculations 
of 

•j 




X" 

id ”*» w o i * r fo 

RC3)= X : utero. 3 in* I *> I of ; Xi 

t— I v 


(see Figure 7, blocks 19,20,21, 22, 23, 24, 25, 26, 27, 28, 29) In order 

to keep reading the data under this form the "read" and "go to" statements 

5 

have been used to read the data from 501 to 10 -500. • (see Figure 7, blocks 
19 and 31) 

(3) Final correction. Read the last 500 data values of the last 
record (i.e. 100th record) and .do the following calculations 

Jr s c- 

£CiJ II x u>j,, ; wten j -frm / -to io 0 ; rr x x c j 

C"t • t-i 


(see Figure 7, blocks 32, 33, 34, 35, 36, 37, 38, 39, 40, and 41) 
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The calculation of the' power spectral density has been divided 
into two parts. First, calculate the raw power spectral density by equa- 
tion (25) (see Figure 7, block 63). Second, calculate the smoothed power 
spectral density by equation (26) (see Figure 7, blocks 66, 67, and 68). 
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FIGURE 7 

Flow Chart of Autocorrelation .Function and Fouer Spectrum 
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CHAPTER IV 

EXAMPLE OF SIMULATED DATA 

For the purpose of' checking the computer result for the autocorrela- 
tion function and power spectrum for the simulated data generated by com- 

t 

puter, the analytical solution for the autocorrelation function and power 
spectrum has to he discussed. From the analytical point of view, we use the 
average value of the simulated data. 

4.1 Autocorrelation and Power Spectrum for a Random Background 

From KORN (1966a) and BENDAT and PIERSOL (1966d), we know that the 
random data are not correlated among themselves. However, there is a corre- 
lation among individual terms. For example, x(t^), x(t^ + t) are uncorrela- 
ted for every At r 0, therefore, autocorrelation function will be a delta 
function at lag time zero, i.e. t = 0 for infinite length of record. The 
power spectrum would be a constant over all the frequency range. This is so 
called white noise. In other words, white noise has a constant power spec- 
tral density, i.e. G(f) = a, R(t) = a5(T> where a is constant. Unfortun- 
ately, such a process for white noise is not physically realizable since the 
variance or R(0) is infinite. This is true for only infinite length of 
record. 

In practice, white noise is approximated by various types of wide- 
band noise, having approximately constant spectral density over a frequency 
band of interest ("band-limited white noise") . From BENDAT and PIERSOL 
(1966e) , bandwidth limited white noise is a random process with a constant 
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power spectrum defined by 


G(f) - a 
= 0 


0 < f Q - (B/2) < f < f + (B/2) 
otherwise 


( 27 ) 


where is the center frequency, and B is the bandwidth. From equation 

(15) it follows that the associated autocorrelation function is 
r di -K/0 , $I»7T8' 

R(t) 


/'* + (%> , $i»tt er , , w „ 

= a cos 2Hf df = aB ( — * — ) 0>s £> c 

J. .» • ’rsr ' 


'f.-ik i 


B 


For the low-pass white noise, f = y , then G(f) becomes 

G(f ) = a 0 < f < B 

= 0 . otherwise 


(28) 


and R(t) = aB ( 


S ;r,2 7Tti? 

2 rre r 


In Figure 8, the frequency information we want is in the low range 
of frequencies. So the low-pass white noise is of interest 1 . The aucocor- 
relation function for low-pass white noise is a sine function which looks 
like the amplitude of the diffraction pattern for the single slit. 
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4.2 Autocorrelation and Power Spectrum Prediction for 
the Average Pulse Signal Above an Average Random Background 

In this section we discuss a square pulse above the random back- 
ground with period T. .If we make this average square pulse signal with a 
100% intensity increase over the background for a time T/2 and off for 
<t:he remaining time T/2, then we ask: What will the autocorrelation func- . 
tion and power spectral density function look like? 

Before analyzing and making this clear, the Autocorrelation theo- 
rem will have to be reexamined. 


Autocorrelation Theorem 

If x(t) has the Fourier transform A(f), then its autocorrelation 
function [ X has the Fourier transform j A (f ) | ^ , i.e. 

J* / Ac-fjJ 2 ^ 2 ^ c df ~ J x*(t)z(t*'irj£fy (16) 

— £*) — od 


This is the unnormalized autocorrelation function with zero mean. 


that 


Suppose there is a rectangle function II (t) which is defined such 

= i i*i<£ 


XOt) 


- o 


hi >y 


The Fourier transformation of the function II (t) is sine (f) = — — — — 

7Tft 

Therefore, using the Autocorrelation theorem, the R(t) is the Fourier trans- 
2 

form of sine f function is the triangle function of unit height and area. 


This function A(t) is defined that 

= I — I ~| ITI< I 

= ■ o m>| 
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These transformations are shown by three^ solid arrows and one 
» , 

dashed arrow in Fig. 9. Any arrow represents a Fourier transformation. 

The lower solid arrow in Fig. 9 indicates that power spectrum is trans- 
formed to the autocorrelation function and the dashed arrow represents the 
inverse transformation. 

From the above transformations we know that a II (t) function for 
the signal will result in a triangular shaped autocorrelation function. 

If you concentrate, the positive side of t axis of II(t), you will have 
a straight line together with two axes which can form 45° right triangle. 

The autocorrelation function of the average square pulse signal above an 
average background will be of a triangular shape along the lag time axis. 

For the purpose of making the autocorrelation function of the 
square pulse signal clear, we show a plot of the average square pulse sig- 
nal above the average background, and its autocorrelation function in 
Fig. 10. Fart (a) of Fig. 10 shows the "square pulse signal above an aver- 
age background. In order to calculate an autocorrelation function for 
zero mean, we use the data in part (b) of Figure 10 which is made by 
shifting the mean value of signal and background. Part (c) of Figure 10 
represents the autocorrelation function for the case of zero mean. It is 
seen that the period of the autocorrelation function is the same as that 
of average signal above the average background. 

Mathematically speaking, the autocorrelation function qualitatively 
can be considered as the product of the same waveform with a lag time t 
between each other. 

(1) When the lag time t equals zero, the normalized autocorrelation 

function must be the maximum value 1 [see Figure 10 (c), a]. 

T 

1(2) When the lag time t = , where T is the period of signal, the 

product is zer.o [see Figure 10 (c) , b] . 
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(a) 



(b) . 



FIGURE 10 
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T 

(3) When the lag time t = y , i.e., out of phase completely, the 
product is minimum value - 1 [see Figure 10 (c), cj . 

^ *y 

Similarly you will have R(~2f )~0 , — -1 , . 

Power Spectrum 

There are two ways to calculate the power spectral density func- 
tion G(f) for the average square pulse above the average background. 

(1) The density function G(f) can be considered as the Fourier 

transform of the autocorrelation function. In other words, the shape of 
2 

G(f) will be sine f function of the positive range when A(t) is the auto- 
correlation function. In this plot there will be several peaks at several 
particular frequencies but the amplitude of these peaks will be gradually 
smaller and then die out for increasing frequency f. 

(2) Ratio of the power of one component to the other: 

From the mathematical derivation equation (13) we have the result 
oo 

g(P) = XW* ~ j where X(t) 2 is referred to as the time averaged 

power of the mean power in X(t). Consider 

Xct; - 

where f = ~ , and a_, a and b are the usual Fourier coefficients, 
n T 0 7 n n 

Then by the orthogonality of the sine and cosine functions 



i 

wWe = C„ = C 


The amplitude, of a certain frequency, f^ can be found out by taking the 

1 2 

area under the peak above the average noise level and equating it to . 
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Theoretically, in the case of a square wave X(t) = ~ 27 
the ratio of the power of the fundamental component to that of first har- 


monic component is 


I 


£c, • K - c * •• c? =(tJ • (t; = ? ■■ i 


Similarly ; £ — 3 $ > 1 etc. 
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4.3 Autocorrelation and Power Spectrum Prediction for the 
Average Pulsar Signal Above an Average Random Background 


The pulsar type signal used is a pulse which is off 35 ms with an 
average counting rate '1 count/ms and on 5 ms with an average counting rate 
2 count/ms. In order to calculate the autocorrelation function with zero 
mean, we have to transform from this [Figure 11 (a)] to the one which forms 
the autocorrelation function with zero mean as shown in (b) of Figure 11. 
Figure 11 (c) indicates the corresponding autocorrelation function for the 
case of (b) . 

For the same reason as discussed for the autocorrelation of an 
average square pulse signal the shape of the autocorrelation function is 
triangular. KORN, (1966b) shows that the minimum value is a negative num- 
ber -ct 2 /N = -1/8 where cr 2 is the variance which is equivalent to autocor- 
relation function at zero lag time, and N is the period divided by the on 
time. ■ For the case (b) Figure 11, 'the triangle has a base of 10 ms. The 
period of autocorrelation is 40 ms which is same as the signal [Figure 11 ' 
(c)] . When the lag time is 40 ms, i.e., the period, the two identical 
waveforms are superimposed completely; i.e., the product of two identical 
data is maximum and this is the case whenever the lag time is a multiple 
number of periods. It is a constant negative number -1/8 for the rest of 
ranges except when r is in the range of 0^5 ms, 35^45 ms, 75^85 ms, etc. 

The simulated data X(t) in Figure 11 (b) has a Fourier series of 


the form: 



od 

+ X! S'* C 


~ ' " 0—1 


T 


) 





FIGURE 
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Using the above result, we can compare the power spectrum calculated by the 
program with the simulated average pulsar signal above an average random 
background . 
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4'. 4 Computer Results for the Autocorrelation Function and the 
Power Spectrum With Several Kinds of Simulated Data 

The analytical solutions have been . examined in the last three sec- 
tions of this chapter. 'In this, section-, I am going to discuss the autocor- 
relation function and power spectrum for the practical data generated by 
the computer. 

.(1) Random Noise 

As I mentioned in Section 4.1, R(x) would be infinite at the origin 
and zero throughout all the range of lag time t in the plot R(x) vs x, and 
G(f) would be constant throughout all the range of frequency f. In other 
words, the autocorrelation function is simply a delta function and the con- 
stant power spectrum is the so-called white noise. 

Practically, all we can have is a definite length of random data. 
Therefore, from equation (28) the autocorrelation function would be an 
extremely large value at the origin, but fluctuates with a small deviation 
about zero throughout the whole range of lag time. [see Figure 12 (a)] The 
power spectrum also will fluctuate with a small deviation about a constant 
mean value throughout the range of low frequency we are interested in. This 
is called low-pass white noise. [see Figure 13 (a)] 

The results of R(x) and G(f) for the simulated random data are shown 
in the Figure 12 (a) and Figure 13 (a) respectively. Note the simulated 
random data has been generated for the average counting rate, 1 count/ms. 

5 

The sample time was 1 ms, and there were 10 data points for total record 
length of 100 sec. 
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(2) Square Pulse Signal Above a Random Background 

The simulated square pulse signal above the random background fluc- 
tuates about the average square pulse above the average random background. 

In other words, simulated data generated by the computer fluctuates about 
the average value which was used in the theoretical analysis in Sections 
4.1, 4.2, and 4,3. For this reason the autocorrelation function at zero 
lag time has an extremely large valu'e which is much larger than that expec- 
ted initially as shown in Figure 11 (c) . 

The data I used is for an average background rate 1 count/ ms for 
20 ms and an increased average rate to 2 counts/ms for 20 ms. In other 
words , this is the pulse 20 ms off and 20 ms on with 100% intensity above 
the random background. In order to calculate the autocorrelation function 
vrith the zero mean, the mean value has to be subtracted from each datum. 

This square wave would have 1 unit amplitude difference between maximum and 
minimum. 

As shown in Section 4.2, the autocorrelation function for this kind 
of data gives a result that is almost the same as shown in Figure 10 (c) 
except at zero lag time. The shape of R(t) is periodic triangular type with 
a period 40 ms which is the same as the period of the square wave. [see 
Figure 12 (b)] 

The power spectrum G(f) of the average square wave data has the lar- 
gest high peak at the fundamental frequency (25 cycles/sec) and a second 

peak at the first harmonic frequency (75 cycles/sec) and a third peak at 

2 

125 cycles/sec', etc. The shape of G(f) is approximately sine f function as 
shown in Figure 13 (b) . 

The power ratio of the fundamental component to the first harmonic 
component is 9/1 and that of the fundamental component to the second harmon- 




FIGURE 12 




50 


ic component is 25/1 theoretically. From the plot of G(f) vs f (Figure 13), 

2 

and Table 3, we calculate the background of G.(f) is 3.008 (counts /HZ). The 
computer results gives the power ratios 8.558/1 and 25.349/1 with the devia- 
tions from the theoretical value of 4.9% and 1.4% respectively. However, 
the value of G(f) beyond the third peak at 175 cycles/sec is not so good 
compared with fhe theoretical value. The method of calculation of the power 
ratios is shown in the Appendix D. 

(3) Pulsar Type Signal Above a Random Background 

These simulated data are the same as for the case of the square 
pulse above a random- background data except that the pulse type signal is 
off 35 ms and on for 5 ms.’ 

Like the case of the square pulse signal above random background, 
the autocorrelation function has an extremely large value at zero lag time. 
This value is different from the one at large: lag times. This extremely 
large value is about ten times larger than the other peak values, because 
two identical fluctuating pulsar waveforms are completely in phase at the 
zero lag time. In the Figure 12 (c>, there are two triangles with peaks at 
40 ms and 80 ms, respectively. The base of the triangle is 10 ms. 

-The calculation of power spectrum can be checked by the power ratio 
method. The power spectrum G(f) of the pulsar type signal above random 
background has the largest peak at the fundamental frequency (25cycle/sec) 
and a second peak at the first harmonic frequency <50 cycle/sec) and a 
third at 75 cycle/sec. [see Figure 13 (c)] 

Theoretically , the power ratio of the fundamental component to the 
second harmonic component is 1.57/1, and that of the first harmonic compo- 
nent to the third harmonic component is 2/1. From Figure 13 and Table 4. 

2 

the power spectrum of the random background of G(f) is 2.248 (counts /HZ). 
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•3. f 6500950 

37 

0. 133 42.42 8 


3.32032967 

39 

0 .115617 090 


3.07000446 

30 

0.1 34091 15 


3.0010824? 

40 

0. 21 716279 


2.° 83433 5 3 

4 l 

0. 27*207 84 3 


2 . -93 59142 

4? 

0.2(1389372 


3.03365154 

43 

0.1 8459249 


3.09661770 

4 4 

0.16136748 


3. 1 1090660 

45 

0. 14315361 


‘3.29 5 387 2 7 

46 

0.11523497 


3 .257 201 19 

47 

0.00401387 


2.990 19814 

4 9 

0.06061661 


2.°?! 43631 

4 0 

0.038434 07 


2.00792942 

50 

0.01875231 


2.92412949 
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' ■- -0.00931841 

52 

-0.02269282 

53 

-0.04548740 

54 

-0.06685668 

•55 

-0.09799665 

'56 

-0. 10921627 

57 

-0. 142442 88 

5 3 

-0.1 5374672 

59 

-0. in] 77997 

60 

-0.2 0400161 ' • 

61 ’ 

-0. 2^-2 6? 869 

6 2 

-0.70773232 

63 

-0.1 87 53499 

6 4 

-0.154 743 03 

65 

-0.1400427? 

6 b 

-0, 1 2049^39 

67 

-0.09132838 

63 ' 

-0.06784606 

69 

-0.04742037 

70 

-0.02699487 

71 

■-0.00150418 

72 

0.01964332 

73 

0.04824098 

74 

0. 06390601 

75 

0.09397078 

76 

'0. 11043715 

77 

0.14000326 

7 8 

0.15429318 

79 

• 0. 17974269 

30 

0.21000032 

• 8 1 

0. 2.322 7 769 

33 

0 . 20423537 

33 

0.1365748? 

8 4 

0.16334553 

8 5 

0.14215624 

86 

0. 11398635 

87 ■ 

0.08794248 

33 

0 . 05977638 

89 

0.03124937 

90 

0.02007224 

91 

-0. 0030S540 

92 

-0.02097471 

' 93 

-0. 0497] 886 

94 

-0.06927502 

o 6 

-0.09534270 

°6 

-0.1124774? ^ 

97 

-0.142158,9 3 

93 

-0.15275878 

99 

-0.187441 17 

100 

-0.20957309 


2.9945 20 3 6 ; 

2.95705427 ! 

2 . 9 ] 1 3845 8 
2,988 74950 
3. 1 1750339 
3.0 76 260 57 
2.96743293 
2.9^A4L445 
2.’) 23310 23 

2 • 5*2 "•>» 1 3 69 * - * 

3.01554394 | 

3.0 30 548 10 ‘ 

2.97146893 j 
2.° 76674 OB | 
3.12944701 * 

3.19570923 \ 

3.10769272 1 

3.06319153 i 
3,05442333 ; 

3.02104282 ; 

2.99955273 i 
3.03165531 ‘ 

3.06587410 
3.08440113 
3.12718863 
3.08937263 ! 

3.02504063 j 
3.04276943 ; 

■2.93475075 j. 
2.63011456 j 
2.91168380 { 

2. 9 806.37 14 j 

2.00754261 { 

3 . Oo'i 570 54 

3.14396381 
3.10764503 ; 

3.03852654 * 

3.01292229 < 

2.96162510 ; 

2.93826694 
3.12558365 
3.13164616 
2.99409485 
3.10326111 : 

3,23718643 
3.05253220 > 

2.99043751 
3.05136583 
2.94631549 
•2.05121822 
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1 

1. 23^77583 


2.24219894 

• -2 

0.08200514 


2.2484064 1 

3 

0.05716614 


2.27999592 

A • 

0.03336975 


2.319 53049 

5 

0.01275801 


3.80827808 

6 

-0.01628971 


5.12577915 

7 

-0.01361207 


’3 .55620395 

9 / 

-0.01756435 


2. 166 447 64 


-0.01518434 


2.28440205 

10 

- 0,00984651 



n. 

- 0.01661853 


4 ,6469998 1 

12 

-0.01158940 


‘3.20044231 

13 

-0.01269277 


2.12567234 

14 

-0.01237036 


2.3 3790234 

15 

-0.012491 75 


3.40310764 

.16 

-0.01921413 


4. 15856171 

17 

-0.01986755 


3.08336353 

L8 

-0.01389924 

* 

2.24613667 

19 

-0.01839143 

\ 

2.30779648 

20 

21 

22 

-('.01372 823 
-0. 01 1 83 049 
-0.02295482 

■ 

3. 10451031 
3.56556892 
2.713518 14 

23 

-0.01185386 

/ 

2. 17718506 

24 

-0.01904688 

2.27599621 

25 

-0 .010061 16 

- ^ 

2.68459225 

.26 

-0.01475369 


2.9 1632175 

27 

-0.01448439 

1 / 

2 .480851 17 

28 

-0.01738782 

* i 

2. 15 56 30 1 1 

29 

-0.01266667 


2.1 /6 39446 

30 

-0.01702931 


’2.39418888 

31 

-0.01823109 


2.5808668 l 

32 

-0.01640008 


2.43372345 

33 

-0.01339928 


2.262019 16 

34 

-0.01488116 


2.25613213 

35 

-0.01738278 


. 2.27603912 

36 

-0.02859438 


2.25998020 

37 

0.01446395 


; 2.2268915 2 

38 

0.02969724 


2.22111607 

39 

0.06513095 


2.21027374 

40 

0.08022875 


2. 18391705 

41 

0. ,1086797.1 


2.20649910 

42 

0.08430517 


2.26772690 

43 

0.05794246 


2.30425739 

44 

. 0.02888179 


2.35.) 138 6 6 

45 

0 .00863192 


2.36390305 

46 

-0.01216547 


2.25334454 

'47 

-0.01331612 


2.1 7950630 

48 

-0.01524717 


2 .260 31399 

49 

-0.01596766 


2.381 1950 7 

'50 

-0.01313644 


2.42072296 
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51 

-0.01592108 

52 

-0.01482943 

53 

-0.01630285 

54 

-0.01557263 

55 

-0.01576539 

56 

-0.02212092 

57 

-0.01728968 

58 

-0.02213512 

59 

-0.01148916 

60 

-0.01 076009 

6 1 

-0.01503152 

62 

-0.01574961 

6 3 

— 0.0 109 U 76 

6 4 

-0.01216391 

6 5 

-0.01090835 

66 

-0.01276767 

67 

-0.01801877 

68 

-0.01 522582 

69 

-0.01611783 

70 ■ . 

-C. 01601792 

71 

-0.01016211 

72 

-0.01476673 

73 

-0.01473686 

74 

-0.01387888 

75 

-0.02071720 

76 

-0.01343445 

77 

0.01087264 

78 

0.03875672 

79 

0.06787103 

80 

0.08175081 

8 1 . 

0.10591346 

82 

0.07978457 

83 

0.05971146 

84 

0.03181953 

85 • 

0.01426323 

86 

-0.015408C9 

87 

-0.01236438 

88 

-0.01426741 

89 

-0.01941834 

90 

-0.01624508 

9 1 

-0.00618494 

92 

-0.01317261 

93 

.-0.01684612 

94 

-0.01159265 

95 - 

-C. 01721993 

96 

-0.01708999 

97 

-0.02087504 

98 

-0.02069632 

99 

-0.0 061 2 LOO" 

100 

-0. 01110787 


I 

! 

i 

i 

<■ 

t 





2. 351.91727 

2.2492CL77 

2.20243549 

2 .27754116 

2.4 1936684 

2.37709427 

2.25894451 

2.26483631 

2.304531 10 

2.4 iuo Uj 26 

2.42852116 

2.32445049 

? .2 59 809 49 

2 .29 706669 

2.38499451 

2 .33301642 

2.25986481 

2.25233269 

2.26955414 

2.31738949 

2.38268757 

2.38652134 

2.32282066 

2 . 280 16090 

2. 20579624 

2.12461281 

2. 14761925 

2.18167877 

2. 17203236 

2. 15672398 

2. 14904022 

2. 23986244 

2.3224229S 

2.28068924 

2.241681 10 

2.18393517 

2.11809349 

2.18995857 

2. 35678959 

2.42445564 

2.31532268 

2.26659775 

2.32386780 

2.37490463 

2 . 3900699 6 

2. 2 50 118 2 6 

2. 13265610 

2. 19956779 

2.41071129 

2.55776882 
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The computer results for the power ratio are 1.47/1 for the fundamental 
component to the second harmonic component, and 1.81/1 for the first har- 
monic component to the third harmonic component. The deviations from the 
theoretical average are 6% for the former case and 9.5% for the latter 
case. 

(4) Smallest Detectable Pulsar Using Autocorrelation 
and Power Spectrum Analysis 

The autocorrelation and power spectrum measurements for the pulsar 
signal with 100% intensity above the random background have been shown in 
section 4.4.3, the result of this analysis shows that the autocorrelation 
function and power spectrum are detectable and predictable. 

If we reduce the intensity of signal, can we still detect the auto- 
correlation function and power spectrum? What is the smallest detectable 
pulsar signal above the random background? Experimentally, for the pulsar 
signal with 25% intensity with respect to the random background, the autocor- 
relation function is not detectable. It is more or less random. (see 
Figure 14) But the power spectrum in the frequency domain is detectable. If 
you plot power spectrum versus frequency in a large scale, you still can see 

c' 

the several peaks at the several expected frequencies (see Figure 15). Of 
course, this curve is not as good as for the case with 100% intensity pulsar 
signal. The values of the power spectrum over all the frequency range are 
not fluctuating very much. It can be imagined that the measurement of auto- 
correlation function and power spectrum will be getting worse and worse for 
reducing the intensity of pulsar signal smaller and smaller. If you keep 
reducing the intensity of signal, finally the result will turn out to be the 
case of random background. Then the autocorrelation function and power spec- 






FIGURE 14 
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•trum would be sine function and approximately a constant with small devia- 
tion, respectively. This result is based upon the short data record of 
length 100 sec and the sampling time of 1 ms. 

It is interesting that a pulsar with 0.5% intensity above random 
background can be detected by use of superimposed’ Epoth Analysis as done 
by LARRY ORWIG for NP 0532. (1971). 
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CHAPTER V 

CONCLUSION AND DISCUSSION 

Autocorrelation and power spectrum technique detects periodic sig- 
nals from the random noise. However, the power spectrum also determines 
the power of each frequency component. For the actual X-ray or gamma ray 
pulsar NP 0532, this method reveals the power of the pulsating component. 

Of course, the power of pulsating component of NP 0532 depends on the 
photon energy range you are interested in. 

The simulated pulsar data for the autocorrelation and power spec- 
trum analysis are done for the situation of the pulsar signal with 100% 
intensit 3 r above the random background. For this data, the result of auto- 
correlation and power spectrum shows the pulsating period or pulsating 
frequency component explicitly. 

Reducing the intensity of the signal, the result of autocorrelation 
and power spectrum analysis is not as good as for the case of a 100% inten- 
sity of pulsar signal. The author has done the case of a 25% intensity of 
signal. The result has been shown in Section 4.4.4. The problem arises 
from shrinking the intensity of signal above the random noise. The reason 
is that the fluctuation of the random background is so large that it buries 
the small intensity signal. We define the fluctuation of the random back- 
ground as the noise. Therefore, how to reduce the noise (fluctuation) is 
our main task. There is one way that can be used to solve this problem. 
That is to increase the number of data points or observation time. The 
distribution of the number of counts within a unit time interval .obeys the 
Poisson distribution. Therefore, the standard deviation a equals square 
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root of the mean value of counts. 

for a given average counting rate, is proportional to square root 
of observation time. However, the signal is proportional to the observa- 
tion time. Then the signal to noise ratio must increase some factor 
because of increasing the observation time. For example, suppose the sig- 


nal with 25% intensity on the random background with average, counting rate 

c 

a and the observation time t, then the signal to noise ratio is — - 

. In the long run, the signal to noise ratio is 


ffr-t 

proportional to the square root of observation time for a given average 
counting rate a. 

The power spectral density calculated could be increased if we re- 
strict the analysis to only the frequency range of interest for the simu- 
lated pulsar in the original data or if we restrict the power spectrum 
analysis to only the lowest frequency components of interest. The total 
power contributed by these low frequency components is equal to the area 
bounded by the relevant portion of' curve of the spectrum and the estimated 
mean level of the- background. The area can be calculated by summing up the 
spectrum multiplied by the frequency resolution for the discrete case i.e. 

p= X . We could reduce the frequency resolution to get a 

/- *• 

higher value of the power spectral density for the frequency of interest. 
The frequency resolution can be expressed as follows: 


2. Ya 

Where m is the maximum lag number and At is the sampling time interval. 

In other words, the frequency resolution can be reduced by increasing the 
sample time interval for a given m. 

Take the case which I have used for an example. A small signal 
above the random background would be able to be detected by increasing the 
sampling time to five times longer than used in the original' analysis. 
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That is, if the sampling time is 5 ms rather than 1 ms the peak in the 
power spectral density plot would just be discernible. For the shorter 
sampling time a 25% signal is just discernible. With a 5 ms sampling time 
a signal of magnitude ^ 1 / would be discernible because the 

power is proportional to the square of the amplitude., 
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APPENDIX A 

FLOW CHART FOR GENERATING SIMULATED RANDOM DATA 


A subroutine RANDU was used to generate the random data, corres- 
ponding to some average rate a. The basic principle and the procedure of 
generating random numbers within each specified time interval At is dis- 
cussed in section 3.1. From systems/360 scientific subroutine package, 
this subroutine is expressed by RANDU' (IX, IY, UNIT). IX is the first 
entry; this must be any odd integer* number with nine or less digits. 

After IX entry, IX should be the previous value of IY computed by this sub- 
routine. IY is a resultant integral random number required for the next 

'■ 31 

entry to this subroutine. The range of this number is between 0 and 2 

UNIT is the resultant uniformly distributed, floating point, random number 

in the range 0 to 1 which is the output of this subroutine. 

The SUBROUTINE RANDU (IX,. IY, UNIT) is shown as follows: 


IY-IX*65539 
IF(IY)5, 6, 6 

5 IY=IY+214 7483 647+1 

6 UNIT = IY 

UNIT = UNIT*.4656613E-9 

RETURN 

END 


By use of the time interval distribution UNIT = e at i or t. =-' 
j x a 

where a is the average rate taken here as 1 count/ms, a sequence of the ran- 
dom, time interval could be generated. stands for the time interval 

between the ith- and the (i-KL)th counts. From this {t^} sequence, number of 
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counts, X can be produced. X is the number of counts each characterized 
m m 

by t. which satisfies the following relation 
3 

w-i < xr.-fc. < m C28> 

Where m represents integers which have to be greater than 1. In other words, 

X^ is equal to j under the equation (28) . 

Physically, t^ is the time interval between two consecutive events. 

The simulated random datum X is the number of events between a particular 

m 

time m and m + At where At = 1 ms. Note the average time interval between 
events (counts) is also 1 ms. 

is nothing but the number of 'Addition 1 . This technique can be 
completed by a ’DO’ loop (see Figure 16, blocks 18 and 19) and a testing 

4 

statement (see Figure 16, block 20). These X data are stored in a nine 

m 

track tape with 800 bits/inch and with the form of variable record length. 
There are 100 records in this tape and each record stores 1000 numbers of 

t 

data. 

For the purpose of checking X^ and t_^, 10th, 20th, 30th ... lOOth's 
record of data are printed out by use of the 'Mod 1 function (see Figure 16 
block 23) and testing statement (see Figure 16, block 24). 



START 

Toi * 

IX=65549 



BEGIN DO LOOP! 
11 1=1,1500 



RANDU 

(IX, IY, UNIT) 


IX=IY 




11 

A \ 

30 LOOPX 



/ . GT . \ . 
/ 30 00 \ • 


BEGIN 
20 1= 

DO LOOP 
1,3000 




14 

X(I) 

l°l 

_N0 P 
0] 

. /DO ] 

A. 15 

.ook 

i 

YES 

,16 

o 



r„ 

BEGIN 1 
22 1= 

30 LOOP 
L,1500 




18 

T=T+TIME (I) 
J=1+IFIX (T) 

22 

19 

X(J)=X(J)+1 


FALSE , n j print j 


LIST=TX 


RETURN TO 
SYSTEM 


/DO Look 

YES 21 


WRITE TO DEV 
25 IN 

INTERNAL FORMAT 
FROM THE LIST ■ 


/ modV— - 

' . NE . Q\ 

~ | FAL SE 
TrINT] 25 


LIST = 
MINX 


_L2L 

PRINT 

25 


l:st=(x(i), 

1=0) 

\29 

PRINT 

26 

30 

,IST=TIME 


200 A 31 


END FILE 
25 


C ™LJ 

RETURN TO 
SYSTEM 



FIGURE .16 

Flow Chart for 
Generating Random 
Background 
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APPENDIX B 

FLOW CHART FOR GENERATING PULSAR SIGNAL ABOVE TIIE RANDOM BACKGROUND 

This is almost the same as in the case of the random background, 
except it calls for using the RANDU SUBROUTINE twice. The random times t 

1 

and Qt_^ were generated for counting rates of lct/ms, and 2cts/ms respec- 
tively. Using the same transformation as before, from Qt_^ to Q^, which 
corresponds to above, gives a signal with intensity of 100% increase 
over the random background. However, a periodic signal above the random 
background can be generated by substituting for the particular values 
of X^. The width of the signal and its period will determine these par- 
ticular numbers of X to be substituted. In Figure 17, blocks 43, 44, 45 

m 

and 46, I = 35, 1000, 40; J = 1, 5; K = K + 1; X(I + J) = Q(K) shows the 

method by which a pulsar signal above background can be produced, i.e. 

this pulsar signal would have an on time of 5 ms and an off-time of 35 ms. 

Similarly, if I = 20, 1000, 40; J = 1, 20 statements correspond to a 

square pulse signal with the same period (40 ms) . This signal has an off- 

time of 20 ms and an on-time of 20 ms. 

In order to test the sum of X in each- record an additional varia- 
nt 

ble M (where M= ) must be considered. This could be printed 

tn-l 

out along with X^ and t^ every 10th record. 




WRITE TO DEV 
25 III 

i INTERNAL FORMAT 
H(OM THE LIST 
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APPENDIX C 
PROGRAMS 

Tho 1*0 nro lUva programs in Uh.1.8 appendix: (I) ilio pcnp.mm of gen- 

erating random background, (2) the program of generating square pulse sig- 
nal with 100% intensity above random background, (3) the program of gener- 
ating pulsar type signal with 100% intensity above random background, 

(4) the program of generating pulsar signal with 25% intensity above ran- 
dom background, and (5) the program of calculating autocorrelation function 
and power spectrum. 
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The Program of Generating Random Background 


FORTRAN IV G LEVEL 13 MAIN. DATE = 


0001 


DIMENSION T IME { 1500 ) ,X(30C0) ,R< 1C5) ,Y(1000) 

0002 


INTEGERS X , Y 

0003 


EQUIVALENCE C-X,Y) 

0004 


IX=65549 

000 5 


DO 200 M I N X = 1 , 100 

0006 


DO 11 1 = 1, 1500 

000 7 


CALL RANDU! IX, I Y, UNI T) 

0008 


IX= IY 

0009 

11 

T IM E ( I )=-( ALUGI UNI T).) 

COLO 


IX =0 . 

00 11 


DO 70 1 = 1, 1500 

0012 

70 

TX=TX+T IME ( I ) 

0013 

* • 

I F( TX .GT .300 0. ) GO TO 50 

0014 


DO 20 1=1, 3000 

0015 

20 

xm=c 

00 lo • 


T=0. 

0017 


DO 22 1=1, 1 50C 

C01S 


T=T+T IME ( I ) 

0019 


J = 1 + IFIX{T) j , 

0020 

22 

X( J ) = X{ J ) + 1 ; 

0021 


WRITE < 2 5 > Y 

0022 


MO D=M I NX- 10* (MINX/ 10) 

0023 


IF(MOD.NE.O) GO TO 20 0 

0024 


PRINT 30., MINX 

0025 

30 

FORMAT { * 1MI NX= * , I 5 ) 

0026 


PRINT 25, (X{I), 1=1, J) 

0027 

25 

FORMAT ( * I*, /{ IX, 101 10) ) 

0028 


PRINT 26, TIME 

0029 

26 

FORMAT! ‘TIME*. /I 1X.1CF10.3) ) 

0030 

200 

CONTINUE 

0031 


END FILE 25 

0032 


STOP 

0033 

50 

PRINT 51, TX 

0034 

51 

FORMAT < 1 TX IS «,I8) 

0035 * 


STOP 

0036 


END 


70285 



The Program of Generating Square Pulse Signal 
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With 100% Intensity Above Random Background 


TRAM 

IV G LEVEL 

19 MAIN DATE 

01 


DIMENSION TIMFU500) ,X(3 000) , R{ 105 ), Y( 1000 ) 

02 


DIMENSION QT I ME ( 1 000} ,Q(2000) 

0 3 


I NTEGER *4 X,Y 

04 


INTEGER* 4 Q 


• 

EQUIVALENCE (X,Y) 

Oft 


I X=65549 

07 


DO 200 MI NX=1 ,100 

08 


DO 11 1=1,1500 • 

09 


CALL RANDU( I X,I Y,UNIT) 

10 


I X=I Y 

11 

11 

TIME ( I )=- (ALOG (UNI T) ) 

12 


TX=0. 

13 


DO 70 1=1,1500 

14 

70 

TX=TX+TI ME ( I ) 

15 


I FI TX.CT.3000. ) GO TO 50 

16 


DO 20 1=1,3000 

17 

. 20 

X 

t-H 

1! 

O 

18 


T=0, 

19 


DO 22 1=1,1500 

20 


T=T+T I ME ( I ) 

21 


J = 1 + IF I X( T) 

22 

22 

X{ J ) = X( J> + 1 

23 


00 1 1=1,1000 

24 


CALL RANDU( I X,I Y,UNIT) 

25 


I X = I Y 

26 

1 

QTIME tn =-(ALOG(UNIT) /2.0) 

? 7 


TX=0.0 

2 3 

• 

on 2 1 = 1,1000 

29 

2 

TX=TX+QTI MF ( I } 

30 


IF { TX.GT.2000 . ) GO TO 50 

31 


on 3 1=1,2000 

32 

3 

0 ( I ) =0 

33 


T=0 . 0 

34 


no 4 1 = 1,1000 • 

35 


T=T+QTTMC (I ) 

36 


J =1 + 1 F I X ( T) • 

37 

4 

Q ( J )=Q( J ) + l 

38 


K = 0 

39 


DO 5 1=20,1000,40 

40 


DO 5 J=1 , 20 

i) 


K=K+ 1 

42 

5 

X(I+J)=0(K) > ‘ '• 

43 


M =0 ' 

4 4 


nn 40 1 = 1,1000 

'+5 

40 

M=M + X ( I ) 

46 


WRITE (25) Y 

47 

• 

MOO = M I NX-10* ( MI NX/10) 

48 


IF(MOD.NE.O) GO TO 200 
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IS 
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TfUN IV G LEVEL 19 


MAIN DATE = 71062 l< 


50 

51 

52 

53 
5 A 

55 

56 

57 

58 

5Q 

60 

61 


PRINT 3 0 t M I N X , M . 

30 FORMAT ( * l.MI NX= ‘,I5 T /,» M=*,I5J 

PRINT 25, < X ( I ) , I =1 ,1000) • 

25 FOR MAT { 'l* ,/(lX,10I10> ) 

PRINT 26, TIME 

26 FORMAT! » TIME * , / (1 X, 1 0F10..3 ) ) 

200 CONTINUE 

END FILE 25 

STOP 

50 PRINT 51, TX 

51 FORMAT ( ' TX IS »,I8) ■ 

STOP 

• END 


The Program of Generating Pulsar Type Signal 
With 100% Intensity Above. Random Background 


73 . 


l 

3RTRAN IV G LEVEL 19 MAIN DATE = 71047 2 

. i 


1001 


DIMENSION TI ME ( 1-500) ,X{3000) ,R<105 )-,YUOOO ) 

iOC 2 


DIMENSION QTIME C300) ,Q{600}' 

5 CO 3 


INTEGER*4 X T Y 

K>04 


INTEGERS Q 

)C 05 

i 

EQUIVALENCE (X,Y> 

'1006 


I X=65549 

.'007 


DO 200 MINX =1,100 

>008 


DO 11 1=1,1500 

>C09 


CALL RANDU< I X,I Y,UNIT) . 

>010 


I X=I Y ■ • 

>011 

11 

TIME ! I ) =-{ALOG(UNIT) ) 

>n 12 • 


TX=0. "* 

>013 


DO 70 1=1,1500 

>0 14 

70 

TX=TX+TIME ( I ) 

K>15 


IF(TX.GT.3000.)G0T0 50 

>016 


DO 20 1=1,3000 • 

>C 17 

2° 

X { I )=0 

>018 


T=0. 

>019 


DO 22 I =1,1500 ■ 

>G2Q 


T=T +TIME (I) 

>0 21 


J = 1 + 1 FI X ( T) 

>022 

22 

X{J)=X(J)+1 

>0 23 


DO 1 1=1,300 

>0 24 


CALL RANDU( I X, I Y,UNI T) 

>0 25 


I X = I Y 

>026 

1 

QTIME (I > =-IvA LOG (UNIT) /2.0) 

>027 


TX=0 ■ 0 

>028 


DO 2 1 = 1,300 • 

3029 

2 

TX=TX+QTIME( I) 

>0 30 


IF (TX.GT. 600.0). GO TO 50- 

>0 31 


DO 3 1=1,600 

>C 32 

3 

Q ( I ) =0 

)0 33 


T=0. 0 

3034 


DO 4 1=1,300- 

30 3 5 


T=T+QTIME( I) 

>0 36 


J=1+IFI X (T) 

30 37 

4 

Q{J)=Q( J>+1 

>038 


K=0 

>039 


DO 5 1=35,1000,40 

3040 


DO 5 J=1 , 5 

>041 


K =K+ 1 

3042 

5 

X ( I+J ) =Q { K) 

>043 


m=6 - • ‘ 

3044 


DO 4Q 1=1,1000 

3C45 

40 

M=M+ X ( I } 

>046 


WRITE (2 5) Y 

>0 47 


MOD =MI NX-10* (MI NX/10) 

3048 


IF ( MOD . NE .0 ) GO TO 200 
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\TRAN IV G LEVEL 10 MAIN DATE = 71047 2 

‘49 PRINT 30, MINX 

50 30 FORMAT <*1MINX= *,I5> 

>51 PRINT 27, M 

>52 27 FORMAT ( *M= *,110) 

>53 PRINT 25, ( X < I ) , 1 = 1 » 100C ) 

554 25 FORMAT! * 1 * , / 1 1 X, 101 10) ) 

555 PRINT 26, TIME 

;,56 26 FORMAT! * TIME 1 ,/tlX,10F10«3.) ) 

251 200 CONTINUE 

> 58 END FILE 25 

059 STOP 

h60 50 PRINT 51 , TX 

b61 • 51 FORMAT ( » TX IS ■ ,18) 

062 • STOP 

p63 END 


The Program of Generating Pulsar Signal 
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With 25 % Intensity Above Random Background 


TRAN IV G LEVEL 19 


MAIN DATE = 71092 


01 

02 

03 

04 

05 

06 
07 
03 

09 

10 

11 ’ 
12 

13 

14 
1 5 
16 
.17 
13 

19 

20 
21 ■ 
22 

23 

24 

25 

26 
27 
23 

29 

30 

31 

32 

33 

34 

35 

36 

37 
33 

39 

40 

41 

42 

43 

44 

45 
)46 
)47 
)48 


DIMENSION TI ME ( 1 500 ) ,X(3000) T R ( 10 5 ) , Y ( 10 00 ) 
, DIMENSION QTIME {10005 ,0(2000) 

INTEGERS X , Y 
INTEGERS Q 
EQUIVALENCE (X,Y) 

I X=65545 

DO 200 M I N X = 1 ,100 

DO 11 1=1,3500 

CALL RANDUII X,I Y,UNIT) 

I X=I Y 

11 TIME (I ) =- ( ALQG (UNI T) ) 

TX=0» 

no 70 1=1,1500 
70 TX=TX+TIMF (I ) 

IFCTX.GT.3000. )G0 TO 50 
no 20 1=1,3000 
20 X( I)=0 
T =0. 

DO 22 1=1,1500 
T=T + T IMF ( I 5 
J =1 + 1 F I X ( T) 

22 X(J)=X(J)+1 
DO 1 1=1,1000 
CALL RANDU(I X,I Y,UNIT) 

.1 X=I Y 

1 QTIME (I) =-( A LOG (UNIT) /l. 25) 

TX=0. 0 

DO 2 1=1*, 1000 

2 TX=TX+QTIME (I ) 

IF ( TX.GT.2000. ) GO TO 50 
DO 3 1=1,2000 
■3 Q ( I ) =0 
T=0. 0 

DO 4 1=1,1000 
• T = T+Q TI ME d ) 

J =1 + 1 F I X { T) 

-"4 Q ( J ) =Q ( J) + 1 
K =0 

DO 5 1=35,1000,40 
DO 5 J=1 , 5 
K =K+ 1 

5 X (I + J 5 =Q ( K) ■ 

M =0 - 

DO 40 1=1,1000 
40 M =M + X ( I ) 

WRITE ( 2 55 Y 

HOD =M I NX- 10* ( MI NX/1 0) , 

I F ( MOD. NE , 0) GO TO 200 
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MAIN DATE = 71092 


49 

90 

51 

52 

53 

54 

5? 

56 

57 

58 

59 

60 
61 


PRINT 30, MINX, M 

30 FORMAT ( * lTfl NX= *,I5,/,* M=',I5) 

•PRINT 25, (X(I) ,1=1 ,1000) 

25 FORMAT! • 1 S / C 1 X , 1 01 1 0) ) 

PRINT 26, TIME 

26 FORMAT ( ‘TIME * , / U X * 1 0F10 .3 ) } 

200 CONTINUE 

END FILE 25 
STOP 

50 PRINT 51, TX 

51 FORMAT ( * TX I S • ,18) 

STOP 

END 


22 / 
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The Program of Calculating Autocorrelation Function 
' and Power Spectrum 



\N IV G LEVEL 10 MAIN DATE * 71095 03/44/ 

DIMENSION X(IOOO) T YU 000) ,A { 1000 ), R { 100 ) ,S ( 100 ), U (100 ), V( 100 ) , 

♦ COSINE (199) , RR (100) 

COMMON /TSAI/X,Y 
INTEGER *4 X ,Y-,A , S , R , AVER 
EQUIVALENCE ( X ( 5 0 1 ) t A ( 1 ) ) 

C CALCULATE INITIAL CORRECTION 

N = 1 0 0 

READ (25) X 

S(1)=0 • / ■ 

R ( 1 ) =0 
DO 1 1=2, N 
S ( I >=0 
R ( I ) =0 
K = I- 1 

. 00 1 J = 1,K 

1 S ( I ) = S ( I ) + X { J ) 

AVER=0 

NUM6ER=0 

DO 2 1=1,500 
A V5R=A VER+ X ( I ) ’ 

DO 2 J=1,N 

2 R(J)=R(J)+X( I)*X(I+J-1) 

NUMBER=NUMBE R^-500 

C MAIN LOOP 

3 READ ( 2 5 ,GND =6) Y 

on 4 1 = 1,1000 

A VER=AVER+A ( I ) 

nn 4 j=i ,n 

4 R { J ) =R ( J > + A ( I ) *A { I + J-l > 

NUMBER =NUMBE R+1000 

00 5 1=1,1000 

5 X ( I ) =Y( I ) 

' GO TO 3 

C FINISH CALCULATION FCR LAST RECORD 

6 DO 7 1=1,500 

7 A VER = AVF.R + A ( I ) 

NUMBER=NUMBE R+500 

A VE =F LOA T (AVER) /F LOAT ( NUMBER) 

A VE2=AVE*AVE 
DO 0 J = 1 , N 
K=501-J ( •*" 

DO 8 1=1, K 

3 R { J ) =R { J ) +A ( I ) ♦A ( I + <J-1 ) 

DO 9 J=2 , N 
K=502-J 

‘ 00 9 I = K-, 500 

9 SU)~S{ J)+A (I )• 

DO 10 J = 1 , N 
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FORTRAN 

00 A 5 
00 A 6 
0047 

00 4 8 
004 9 
00^0 
0051 
00 52 
00 5 3 
00 54 
0055 
00 56 
00 57 
00 5 8 

0059 

0060 
0061 
0062 

0063 

0064 

00 6 5 
0066 


IV G LEVEL 19 MAIN DATE = .71095 

10 RR{ J}-(R ( J) + S( J) *AVE-AV£2* (NUMBER+J-1 1)/ (NUMBER-J + 1) 
PRINT 51 .NUMBER, AVER, AVE2 ,R ,S 
51 FORMAT <2X,2I 15,F20.10,/ (2X,8I15 ) ) 

C CALCULATE SPECTRAL DENSITY 

P I =3. 141 59265358579 
DO 12 1=0,198 

. 12 COSINF (1+1) =CGSU*PI/99.0) 

DO 13 K=1 ,N ' 

' KL = MOD ( K , 2 ) 

KM=MOD ( K+ 1 ,2) 

V ( K ) -2* 0# ( RR (1)+KL*RR(100) -KM#RR( 100 ) ) 

DO 13 L=2 ,99 
M = { L- 1 ) * ( K- 1 ) 

IF(M.GT.198) M=MOD ( M , 1 98 ) 

13 V(K ) =V (K )+RR { L) *CGSI NE { M+l ) *4.0 
U(1)=0.5*(V(1)+V(2) ) 

DO 14 I = 2 , N 

14 U( I )=C.25*V( I-D+O. 5*V(n+0.25*V ( 1 + 1) 

U(100)=0. 5*( V(99) + V(100) ) 

PRINT -50, ( {I » R R { I ) ,V(I) ,U(I) ,1 = 1 , N ) ) 

50 FORMAT t • 1 1 , / »4X, »I • ,9X, ' R' ,11X, » V' , 11X »' U» ,//,/( 2X , 15, 3 
* ( 2X,F 20. 8) > ) 

STOP 

END 



79 


APPENDIX D 

CALCULATIONS OF THE POWER RATIOS 

The power contributed by each frequency component is proportional 

to the square of its ■ amplitude. It may be calculated from the area under 

the .curve of the relevant portion of the power spectrum i.e. that area 

bounded by the curve of the spectrum between specified frequency limits 

and the estimated mean level of the noise background. 

The large value peaks were analyzed in this fashion. Each peak 

was formed from three consecutive frequencies. This was because the data 

had been obtained in a digital form. 

In this experiment the estimated magnitude of the power of the noise 

in each frequency component was obtained from the mean power for the fre- 

2 

quencies 250 'v 495 HZ. This was calculated to be 3.008 (counts /HZ) for 
Table 3 and 2.248 (counts 2 /HZ) from Table 4. 

For square pulse signals above a random background, the power ratio 
of the fundamental component to the first harmonic component is 9/1 and 
that of the fundamental component to the second harmonic component is 25/1. 
These values were calculated as 8.558/1 and 25.349/1. 

The pulsar analysis gives power ratios of the fundamental frequency 
to the second harmonic component to be 1.57/1 and that of the first" compo- 
nent to the third'as 2/1. The calculated obtained values were 1.47/1 and 


1.81/1. 



